Introduction à la modélisation statistique bayésienne

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 & 4: Modèle de régression linéaire

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Modèles multi-niveaux

Cours 9: Introduction à brms

Cours 10: Data Hackaton

On s'intéresse aux différences de taille entre hommes et femmes. On va mesurer 100 femmes et 100 hommes.

men <- rnorm(100, 175, 10) # 100 tailles d'hommes
women <- rnorm(100, 170, 10) # 100 tailles de femmes

t.test(men, women)

    Welch Two Sample t-test

data:  men and women
t = 2.1475, df = 197.71, p-value = 0.03297
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2655634 6.2364837
sample estimates:
mean of x mean of y 
 173.0157  169.7646 

On va simuler des t-valeurs issues de données générées selon l'hypothèse de non-différence entre hommes et femmes.

nSims <- 1e4 # number of simulations
t <- rep(NA, nSims) # initialising an empty vector

for (i in 1:nSims) {

    men2 <- rnorm(100, 170, 10)
    women2 <- rnorm(100, 170, 10)
    t[i] <- t.test(men2, women2)$statistic

}

t <- replicate(nSims, t.test(rnorm(100, 170, 10), rnorm(100, 170, 10) )$statistic)

data.frame(t = t) %>%
    ggplot(aes(x = t) ) +
    geom_histogram() +
    theme_bw(base_size = 20)

data.frame(x = c(-5, 5) ) %>%
    ggplot(aes(x = x) ) +
    stat_function(fun = dt, args = list(df = t.test(men, women)$parameter), size = 1.5) +
    theme_bw(base_size = 20)

abs(qt(0.05 / 2, df = t.test(men, women)$parameter) ) # two-sided critical t-value

[1] 1.972035

alpha <- .05
abs(qt(alpha / 2, df = t.test(men, women)$parameter) ) # two-sided critical t-value

[1] 1.972035

tobs <- t.test(men, women)$statistic # observed t-value
tobs %>% as.numeric

[1] 2.147452

A p-value is simply a tail area (an integral), under the distribution of test statistics under the null hypothesis. It gives the probability of observing the data we observed (or more extreme data), given that the null hypothesis is true.

\[ p[\mathbf{t}(\mathbf{x}^{\text{rep}}|H_{0}) \geq t(x)] \]

t.test(men, women)$p.value

[1] 0.03297456

tvalue <- abs(t.test(men, women)$statistic)
df <- t.test(men, women)$parameter

2 * integrate(dt, tvalue, Inf, df = df)$value

[1] 0.03297456

2 * (1 - pt(abs(t.test(men, women)$statistic), t.test(men, women)$parameter) )

         t 
0.03297456 

Comparaison de deux modèles:

• \( H_{0}: \mu_{1} = \mu_{2} \rightarrow \delta = 0 \)
• \( H_{1}: \mu_{1} \neq \mu_{2} \rightarrow \delta \neq 0 \)

\[ \underbrace{\dfrac{p(H_{0}|D)}{p(H_{1}|D)}}_{posterior\ odds} = \underbrace{\dfrac{p(D|H_{0})}{p(D|H_{1})}}_{Bayes\ factor} \times \underbrace{\dfrac{p(H_{0})}{p(H_{1})}}_{prior\ odds} \]

\[ \text{evidence}\ = p(D|H) = \int p(\theta|H) p(D|\theta,H) \text{d}\theta \]

L'évidence en faveur d'un modèle correspond à la probabilité marginale d'un modèle (le dénominateur du théorème de Bayes), c'est à dire à la likelihood moyennée sur le prior… Ce qui fait du Bayes Factor un ratio de likelihood pondéré (par le prior).

On lance une pièce 100 fois et on essaye d'estimer la probabilité \( \theta \) (le biais de la pièce) d'obtenir face. On compare deux modèles qui diffèrent par leur prior sur \( \theta \).

\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]

\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]

\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]

\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]

\[ \text{BF}_{12} = \dfrac{p(D|H_{1})}{p(D|H_{2})} = \dfrac{\int p(\theta|H_{1}) p(D|\theta,H_{1}) \text{d}\theta}{\int p(\theta|H_{2}) p(D|\theta,H_{2}) \text{d}\theta} = \dfrac{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(6, 10)\text{d}\theta}{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(20, 12) \text{d}\theta} \]

Attention à ne pas interpréter le BF comme un posterior odds…

Le BF est un facteur qui nous indique de combien nos prior odds doivent changer, au vue des données. Il ne nous dit pas quelle est l'hypothèse la plus probable, sachant les données (sauf si les prior odds sont 1:1).

• \( H_{0} \): there is no such thing as precognition
• \( H_{1} \): precognition exists

On fait une expérience et on calcule un \( BF_{10} = 27 \). Quelle est la probabilité a posteriori de H1 ?

\[ \underbrace{\dfrac{p(H_{1}|D)}{p(H_{0}|D)}}_{posterior\ odds} = \underbrace{\dfrac{27}{1}}_{Bayes\ factor} \times \underbrace{\dfrac{1}{1000}}_{prior\ odds} = \dfrac{27}{1000} = 0.027 \]

library(BayesFactor)
ttestBF(men, women)

Bayes factor analysis
--------------
[1] Alt., r=0.707 : 1.312345 ±0%

Against denominator:
  Null, mu1-mu2 = 0 
---
Bayes factor type: BFindepSample, JZS

Prior JZS, en référence à Jeffreys (1961), et Zellner & Siow (1980).

data.frame(x = c(-10, 10) ) %>%
    ggplot(aes(x = x) ) +
    stat_function(
        fun = dcauchy,
        args = list(location = 0, scale = sqrt(2) / 2), size = 1.5) +
    theme_bw(base_size = 20)

\[ \mathrm{Cauchy} = \dfrac{1}{\pi \gamma \Big[ 1 + \big(\frac{x-x_{0}}{\gamma}\big)^{2}\Big]} \]

Dans le cadre bayésien, le terme de prior peut faire référence soit:

• à la probabilité a priori d'un modèle (par rapport à un autre modèle), c'est à dire \( p(M_{i}) \). Voir ce blogpost.

• ou aux prédictions du modèle, e.g., \( \beta \sim \mathrm{Normal}(0,10) \). Voir ce blogpost.

Deux problèmes récurrents à éviter en modélisation: l'overffitting et l'underfitting. Comment s'en sortir ?

• Utiliser des priors pour régulariser, pour contraindre le posterior (i.e., accorder moins de poids à la likelihood)
• Utiliser des critères d'information (e.g., AIC, WAIC)

\[ R^{2} = \frac{\text{var}(\text{outcome}) - \text{var}(\text{residuals})}{\text{var}(\text{outcome})} = 1 - \frac{\text{var}(\text{residuals})}{\text{var}(\text{outcome})} \]

ppnames <- c("afarensis","africanus","habilis","boisei",
        "rudolfensis","ergaster","sapiens")
brainvolcc <- c(438, 452, 612, 521, 752, 871, 1350)
masskg <- c(37.0, 35.5, 34.5, 41.5, 55.5, 61.0, 53.5)

d <- data.frame(species = ppnames, brain = brainvolcc, mass = masskg)

d %>%
    ggplot(aes(x = mass, y = brain, label = species) ) +
    geom_point() +
    ggrepel::geom_label_repel(hjust = 0, nudge_y = 50, size = 5) +
    theme_bw(base_size = 20) + xlim(30, 70)

mod1.1 <- lm(brain ~ mass, data = d)
(var(d$brain) - var(residuals(mod1.1) ) ) / var(d$brain)

[1] 0.490158

mod1.2 <- lm(brain ~ mass + I(mass^2), data = d)
(var(d$brain) - var(residuals(mod1.2) ) ) / var(d$brain)

[1] 0.5359967

mod1.3 <- lm(brain ~ mass + I(mass^2) + I(mass^3), data = d)
(var(d$brain) - var(residuals(mod1.3) ) ) / var(d$brain)

[1] 0.6797736

\[ \begin{align} v_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha \\ \end{align} \]

mod1.7 <- lm(brain ~ 1, data = d)

On voudrait mesurer la distance entre notre modèle et le full model (i.e., la “réalité”, la nature)…

De combien notre incertitude est-elle réduite en apprenant un outcome supplémentaire ? Cette réduction est la définition de l'information.

Quelle mesure pour l'incertitude ? S'il existe \( n \) événements possibles, et que chaque événement \( i \) a pour probabilité \( p_{i} \), alors une mesure de l'incertitude est:

\[ H(p) = - \text{E log}(p_{i}) = - \sum_{i=1}^{n}p_{i} \text{log}(p_{i}) \]

En d'autres termes, l'incertitude contenue dans une distribution de probabilités est la log-probabilité moyenne d'un événement.

Exemple de prédiction météorologique. Si on imagine que la probabilité qu'il pleuve ou qu'il fasse beau (à Grenoble) sont \( p_{1}=0.3 \) et \( p_{2}=0.7 \).

Alors, \( H(p) = - (p_{1} \text{log}(p_{1}) + p_{2} \text{log}(p_{2})) \approx 0.61 \).

p <- c(0.3, 0.7)
- sum(p * log(p) )

[1] 0.6108643

Imaginons que nous habitions à Abu Dabi, et que la probabilité qu'il pleuve ou qu'il fasse beau sont \( p_{1}=0.01 \) et \( p_{2}=0.99 \).

p <- c(0.01, 0.99)
- sum(p * log(p) )

[1] 0.05600153

On a donc un moyen de quantifier l'incertitude. Comment utiliser cette mesure pour quantifier la distance entre notre modèle et la réalité ?

Divergence: incertitude ajoutée par l'utilisation d'une distribution de probabilités pour décrire… une autre distribution de probabilités (Kullback-Leibler divergence, ou entropie relative).

\[ D_{KL}(p,q) = \sum_{i} p_{i}\big(\text{log}(p_{i}) - \text{log}(q_{i})\big) = \sum_{i} p_{i} \text{log}\bigg(\frac{p_{i}}{q_{i}}\bigg) \]

\[ D_{KL}(p,q) = \sum_{i} p_{i}\big(\text{log}(p_{i}) - \text{log}(q_{i})\big) = \sum_{i} p_{i} \text{log}\bigg(\frac{p_{i}}{q_{i}}\bigg) \]

Par exemple, supposons que la véritable distribution de nos événements (soleil vs pluie) soit \( p_{1}=0.3 \) et \( p_{2}=0.7 \). Si nous pensons plutôt que ces événements arrivent avec une probabilité \( q_{1}=0.25 \) et \( q_{2}=0.75 \), quelle quantité d'incertitude avons-nous ajoutée ?

p <- c(0.3, 0.7)
q <- c(0.25, 0.75)

sum(p * log(p / q) )

[1] 0.006401457

sum(q * log(q / p) )

[1] 0.006164264

Entropie croisée: \( H(p,q) = \sum_{i} p_{i} \log (q_{i}) \)

sum(p * (log(q) ) )

La Divergence est définie comme l'entropie additionnelle ajoutée en utilisant \( q \) pour décrire \( p \).

\[ \begin{align} D_{KL}(p,q) &= H(p,q) - H(p) \\ &= - \sum_{i} p_{i} \log(q_{i}) - \big( - \sum_{i} p_{i} \log(p_{i}) \big) \\ &= - \sum_{i} p_{i} \big(\log(q_{i}) - \log(p_{i}) \big) \\ \end{align} \]

- sum (p * (log(q) - log(p) ) )

[1] 0.006401457

OK, mais nous ne connaissons pas la distribution target (la réalité), à quoi cela peut donc nous servir ?

Astuce: si nous comparons deux modèles, \( q \) et \( r \), pour approximer \( p \), nous allons comparer leurs divergences… Et donc \( \text{E} \ \text{log}(p_{i}) \) sera la même quantité pour les deux modèles !

On peut donc utiliser \( \text{E} \ \text{log}(q_{i}) \) et \( \text{E} \ \text{log}(r_{i}) \) comme estimateurs de la distance entre chaque modèle et notre distribution cible. On a donc seulement besoin de la log-probabilité moyenne des modèles. Comme on ne connaît pas la distribution cible, cela veut dire qu'on ne peut pas interpréter cette quantité en termes absolus mais seulement en termes relatifs. Ce qui nous intéresse c'est \( \text{E} \ \text{log}(q_{i}) - \text{E} \ \text{log}(r_{i}) \).

Pour approximer la valeur de \( \text{E} \ \log(q_{i}) \), on peut utiliser la déviance d'un modèle, qui est une mesure du fit relatif du modèle.

\[ D(q) = -2 \sum_{i} \log(q_{i}) \]

où \( i \) indice chaque observation et \( q_{i} \) est la likelihood de chaque observation.

d$mass.s <- scale(d$mass)
mod1.8 <- lm(brain ~ mass.s, data = d)

-2 * logLik(mod1.8) # compute deviance

'log Lik.' 94.92499 (df=3)

# extracting model's coefficients

alpha <- coef(mod1.8)[1]
beta <- coef(mod1.8)[2]

# computing the log-likelihood

ll <- sum(dnorm(
    d$brain,
    mean = alpha + beta * d$mass.,
    sd = sd(residuals(mod1.8) ),
    log = TRUE)
    )

# computing the deviance

(-2) * ll

[1] 95.00404

La déviance a le même problème que le \( R^{2} \), lorsqu'elle est calculée sur l'échantillon observé. Dans ce cas, on l'appelle déviance in-sample.

Si on est intéressé par les capacités de prédiction de notre modèle, nous pouvons calculer la déviance du modèle sur de nouvelles données… qu'on appellera dans ce cas déviance out-of-sample. Cela revient à se demander si notre modèle est performant pour prédire de nouvelles données.

Imaginons que nous disposions d'un échantillon de taille \( N \), que nous appellerons échantillon d'apprentissage (training). Nous pouvons calculer la déviance du modèle sur cet échantillon (\( D_{train} \) ou \( D_{in} \)). Si nous acquérons ensuite un nouvel échantillon de taille \( N \) issu du même processus de génération de données (que nous appellerons échantillon de test), nous pouvons calculer une déviance sur ce nouvel échantillon, en utilisant les paramètres estimés avec l'échantillon d'entraînement (que nous appellerons \( D_{test} \) ou \( D_{out} \)).

\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i},1) \\ \mu_{i} &= (0.15) x_{1,i} - (0.4) x_{2,i} \\ \end{align} \]

On a réalisé ce processus 10.000 fois pour cinq modèles de régression linéaire de complexité croissante. Les points bleus représentent la déviance calculée sur l'échantillon d'apprentissage et les points noirs la déviance calculée sur l'échantillon de test.

Une autre manière de lutter contre l'overfitting est d'utiliser des priors sceptiques qui vont venir ralentir l'apprentissage réalisé sur les données (i.e., accorder plus de poids au prior).

Comment décider de la précision du prior ? Est-ce que le prior est “assez” régularisateur ou pas ?

On peut diviser le jeu de données en deux parties (training et test) afin de choisir le prior qui produit la déviance out-of-sample la plus faible. On appelle cette stratégie la cross-validation.

On mesure ici la différence entre la déviance in-sample (en bleu) et la déviance out-of-sample (en noir). On remarque que la déviance out-of-sample est presque exactement égale à la déviance in-sample, plus deux fois le nombre de paramètres du modèle…

L'AIC fournit une approximation de la déviance out-of-sample:

\[ \text{AIC} = D_{train} + 2p \]

où \( p \) est le nombre de paramètres libres (i.e., à estimer) dans le modèle. L'AIC donne donc une approximation des capacités de prédiction du modèle.

NB: l'AIC fonctionne bien uniquement quand le nombre d'observations \( N \) est largement supérieur au nombre de paramètres \( p \). Dans le cas contraire, on utilise plutôt l'AICc (corrected AIC, voir Burnham & Anderson, 2002; 2004).

Une autre condition de l'AIC est que les priors soient plats (i.e., peu informatifs) ou dépassés par la likelihood (on a beaucoup de données). Le DIC est un indice qui ne requiert pas cette condition, en s'accommodant de priors informatifs.

Le DIC est calculé à partir de la distribution a posteriori de la déviance calculée sur l'échantillon d'apprentissage (i.e., \( D_{train} \)).

\[ \text{DIC} = \bar{D} + (\bar{D} - \widehat{D} ) = \widehat{D} + p_{D} \]

où \( \bar{D} \) est la moyenne de la distribution a posteriori \( D \) calculée pour chaque valeur de paramètre échantillonnée, et \( \widehat{D} \) la déviance calculée à la moyenne de la distribution a posteriori. La différence \( \bar{D} - \widehat{D} = p_{D} \) est analogue au nombre de paramètres utilisé dans le calcul de l'AIC.

La fonction rethinking::DIC permet de le calculer.

Une condition d'application de l'AIC et du DIC est que la distribution a posteriori soit une distribution gaussienne multivariée. Le WAIC relache cette condition, en étant souvent plus précis que le DIC.

Un aspect important du WAIC est qu'il est dit pointwise, c'est à dire qu'il considère l'imprécision de prédiction point par point (donnée par donnée), indépendemment pour chaque observation.

On va commencer par calculer la log-pointwise-predictive-density (lppd), définie par:

\[ \text{lppd} = \sum_{i=1}^{N} \log \Pr(y_{i}) \]

où \( \Pr(y_{i}) \) représente la likelihood moyenne de l'observation \( i \) dans le training sample.

En français: la log-densité prédictive point par point est la somme du log de la likelihood moyenne de chaque observation. Il s'agit de l'analogue point par point de la déviance, moyenné sur toute la distribution a posteriori.

library(rethinking)
data(cars)

m <- rethinking::map(alist(
        dist ~ dnorm(mu, sigma),
        mu <- a + b * speed,
        a ~ dnorm(0, 100),
        b ~ dnorm(0, 10),
        sigma ~ dunif(0, 30) ),
        data = cars
    )

post <- extract.samples(m, n = 1e3)

head(post)

          a        b    sigma
1 -20.74432 4.068218 16.08980
2 -20.34022 4.146951 12.95898
3 -17.60990 3.990101 17.70716
4 -17.00744 3.963182 17.01450
5 -26.66950 4.280532 16.19768
6 -19.28354 3.969091 15.49696

\[ \text{lppd} = \sum_{i=1}^{N} \log \Pr(y_{i}) \]

ll <- sapply(
    1:nrow(post),
    function(s) {

        mu <- post$a[s] + post$b[s] * cars$speed
        dnorm(cars$dist, mu, post$sigma[s], log = TRUE)

        }
    )

lppd <- sapply(1:nrow(cars), function(i) log_sum_exp(ll[i, ]) - log(nrow(post) ) )

sum(lppd)

[1] -206.5098

La deuxième partie du calcul du WAIC est le nombre de paramètres effectif, \( p_{WAIC} \). On définit \( V(y_{i}) \) comme la variance de la log-likelihood pour chaque observation \( i \) de l'échantillon d'entraînement.

\[ p_{\text{WAIC}} = \sum_{i=1}^{N} V(y_{i}) \]

pWAIC <- sapply(1:nrow(cars), function(i) var(ll[i, ]) )

Ensuite, le WAIC est défini par:

\[ \text{WAIC} = -2(\text{lppd} - p_{\text{WAIC}}) \]

(WAIC <- (-2) * (sum(lppd) - sum(pWAIC) ) )

[1] 421.0726

Le WAIC est aussi un estimateur de la déviance out-of-sample. La fonction WAIC du package rethinking permet de le calculer pour des modèles construits avec map ou map2stan.

Le DIC et le WAIC peuvent être conceptualisés (au même titre que l'AIC) comme des approximations de la déviance out-of-sample. On remarque que le WAIC produit des approximations plus précises que le DIC, et que l'usage combiné des critères d'information et des priors régularisateurs permet une meilleure approximation de la déviance out-of-sample.

• Sélection de modèle. On choisit le meilleur modèle un utilisant un des outils présentés et on base nos conclusions sur les paramètres estimés par ce meilleur modèle.

• Comparaison de modèles. On utilise DIC/WAIC mais aussi les outils de posterior predictive checking discutés précédemment, pour chaque modèle.

• Moyennage de modèles. On va construire des posterior predictive checks qui exploitent ce qu'on sait des capacités de prédiction de chaque modèle (via le DIC/WAIC).

Nous allons essayer de prédire les kg par gramme de lait (kcal.per.g) avec les prédicteurs neocortex et le logarithme de mass. Nous allons ensuite fitter 4 modèles qui correspondent aux 4 combinaisons possibles de prédicteurs et les comparer en utilisant le WAIC.

data(milk)

d <- milk[complete.cases(milk), ] # removing NAs
d$neocortex <- d$neocortex.perc / 100 # rescaling explanatory variable

head(d)

              clade            species kcal.per.g perc.fat perc.protein
1     Strepsirrhine     Eulemur fulvus       0.49    16.60        15.42
6  New World Monkey Alouatta seniculus       0.47    21.22        23.58
7  New World Monkey         A palliata       0.56    29.66        23.46
8  New World Monkey       Cebus apella       0.89    53.41        15.80
10 New World Monkey         S sciureus       0.92    50.58        22.33
11 New World Monkey   Cebuella pygmaea       0.80    41.35        20.85
   perc.lactose mass neocortex.perc neocortex
1         67.98 1.95          55.16    0.5516
6         55.20 5.25          64.54    0.6454
7         46.88 5.37          64.54    0.6454
8         30.79 2.51          67.64    0.6764
10        27.09 0.68          68.85    0.6885
11        37.80 0.12          58.85    0.5885

a.start <- mean(d$kcal.per.g)
sigma.start <- log(sd(d$kcal.per.g) )

mod2.1 <- rethinking::map(alist(
        kcal.per.g ~ dnorm(a, exp(log.sigma) ) ),
        data = d,
        start = list(a = a.start, log.sigma = sigma.start) )

mod2.2 <- rethinking::map(alist(
        kcal.per.g ~ dnorm(mu, exp(log.sigma) ),
        mu <- a + bn * neocortex),
        data = d,
        start = list(a = a.start, bn = 0, log.sigma = sigma.start) )

mod2.3 <- rethinking::map(alist(
        kcal.per.g ~ dnorm(mu, exp(log.sigma) ),
        mu <- a + bm * log(mass) ),
        data = d,
        start = list(a = a.start, bm = 0, log.sigma = sigma.start) )

mod2.4 <- rethinking::map(alist(
        kcal.per.g ~ dnorm(mu, exp(log.sigma) ),
        mu <- a + bn * neocortex + bm * log(mass) ),
        data = d,
        start = list(a = a.start, bn = 0, bm = 0, log.sigma = sigma.start) )

(milk.models <- compare(mod2.1, mod2.2, mod2.3, mod2.4) )

        WAIC pWAIC dWAIC weight   SE  dSE
mod2.4 -14.7   5.0   0.0   0.92 7.91   NA
mod2.1  -8.4   1.8   6.3   0.04 4.52 7.41
mod2.3  -7.8   3.0   6.8   0.03 5.67 5.35
mod2.2  -6.2   2.9   8.4   0.01 4.28 7.79

plot(milk.models, SE = TRUE, dSE = TRUE)

Le poids d'un modèle est une estimation de la probabilité que ce modèle fera les meilleures prédictions possibles sur un nouveau jeu de données, conditionnellement au set de modèles considéré.

\[ w_{i} = \frac{\exp (-\frac{1}{2}\text{dWAIC}_{i})}{\sum_{j=1}^{m} \exp (-\frac{1}{2}\text{dWAIC}_{j})} \]

Cette formule “remet” simplement le WAIC à l'échelle des probabilités. Le modèle mod2.1 a un poids de \( 0.92 \), qui le place en tête du jeu de modèles. N'oublions cependant pas que nous disposons seulement de 12 observations…

coeftab(mod2.1, mod2.2, mod2.3, mod2.4)

          mod2.1  mod2.2  mod2.3  mod2.4 
a            0.66    0.35    0.71   -1.09
log.sigma   -1.79   -1.80   -1.85   -2.16
bn             NA    0.45      NA    2.79
bm             NA      NA   -0.03   -0.10
nobs           17      17      17      17

plot(coeftab(mod2.1, mod2.2, mod2.3, mod2.4) )

Pourquoi ne conserver uniquement le premier modèle et oublier les autres ? Une autre stratégie consisterait à pondérer les prédictions des modèles par leurs poids respectifs. C'est ce qu'on appelle le moyennage de modèles (model averaging).

• Calculer le WAIC de chaque modèle
• Calculer le poids de chaque modèle
• Simuler des données à partir de chaque modèle
• Combiner ces valeurs simulées dans un ensemble de prédictions, en utilisant les poids des modèles comme proportions

nc.seq <- seq(from = 0.5, to = 0.8, length.out = 30)

d.predict <- list(
    kcal.per.g = rep(0, 30), # empty outcome
    neocortex = nc.seq,      # sequence of neocortex
    mass = rep(4.5, 30)      # average mass
    )

milk.ensemble <- ensemble(mod2.1, mod2.2, mod2.3, mod2.4, data = d.predict, refresh = 0)
mu <- apply(milk.ensemble$link, 2, mean)
mu.HPDI <- t(apply(milk.ensemble$link, 2, HPDI) )

Voici les prédictions de tous les modèles considérés, pondérés par leur poids respectif. Comme le modèle mod2.4 concentrait quasiment tout le poids, il fait sens que cette prédiction moyennée soit similaire aux prédictions du modèle mod2.4.

data(Howell1)

d <-
    Howell1 %>%
    mutate(age = scale(age) )

set.seed(1000)
i <- sample(1:nrow(d), size = nrow(d) / 2)

d1 <- d[i, ] # training sample
d2 <- d[-i, ] # test sample

Nous avons maintenant deux dataframes, de 272 lignes chacune. On va utiliser d1 pour fitter nos modèles et d2 pour les évaluer.

Soit \( h_{i} \) les valeurs de taille et \( x_{i} \) les valeurs centrées d'âge, sur la ligne \( i \). Construisez les modèles suivants avec d1, en utilisant rethinking::map et des priors faiblement régularisateurs.

\[ \begin{align} \mathcal{M_{1}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} \\ \mathcal{M_{2}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} \\ \mathcal{M_{3}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} \\ \mathcal{M_{4}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} x_{i}^{4} \\ \mathcal{M_{5}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} x_{i}^{4} + \beta_{5} x_{i}^{5} \\ \mathcal{M_{6}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} x_{i}^{4} + \beta_{5} x_{i}^{5} + \beta_{6} x_{i}^{6} \\ \end{align} \]

1. Comparer ces modèles en utilisant le WAIC. Comparer les rangs des modèles et leurs poids.
2. Pour chaque modèle, produire un plot de la moyenne estimée et l'intervalle de confiance à 97% de la moyenne, surimposée aux données brutes. Comment ces prédictions diffèrent-elles selon les modèles ?
3. Faire un plot des prédictions moyennées sur tous les modèles (sur les trois meilleurs). En quoi ces prédictions diffèrent-elles des prédictions du modèle avec le plus petit WAIC ?
4. Calculer la déviance out-of-sample pour chaque modèle. Vous pouvez calculer la log-likelihood des tailles en utilisant sum(dnorm(d2$height, mu, sigma, log = TRUE)) où mu est un vecteur de moyennes prédites, et sigma est l'écart type produit par map.
5. Comparer les déviances obtenues à la question précédente aux valeurs de WAIC. Basé sur les déviances obtenues, quel modèle fait les meilleures prédictions ? Est-ce que le WAIC est un bon estimateur la déviance out-of-sample ?

alpha.start <- mean(d1$height)
sigma.start <- sd(d1$height)

mod3.1 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age,
        c(alpha, beta.1) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

mod3.2 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2,
        c(alpha, beta.1, beta.2) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

mod3.3 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2 + beta.3*age^3,
        c(alpha, beta.1, beta.2, beta.3) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

mod3.4 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2 + beta.3*age^3 + beta.4*age^4,
        c(alpha, beta.1, beta.2, beta.3, beta.4) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

mod3.5 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2 + beta.3*age^3 + 
                beta.4*age^4 + beta.5*age^5,
        c(alpha, beta.1, beta.2, beta.3, beta.4, beta.5) ~ dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

mod3.6 <- rethinking::map(alist(
        height ~ dnorm(mean = mu, sd = sigma),
        mu <- alpha + beta.1*age + beta.2*age^2 + beta.3*age^3 +
                beta.4*age^4 +
                beta.5*age^5 + beta.6*age^6,
        c(alpha, beta.1, beta.2, beta.3, beta.4, beta.5, beta.6) ~dnorm(mean = 0, sd = 100),
        sigma ~ dunif(min = 0, max = 50) ),
        data = d1,
        start = list(alpha = alpha.start) )

compare(mod3.1, mod3.2, mod3.3, mod3.4, mod3.5, mod3.6)

         WAIC pWAIC dWAIC weight    SE   dSE
mod3.4 1926.0   5.6   0.0   0.57 25.44    NA
mod3.5 1927.5   6.4   1.4   0.28 25.32  1.15
mod3.6 1928.6   7.8   2.6   0.16 25.04  2.89
mod3.3 1952.3   5.4  26.2   0.00 24.19 10.85
mod3.2 2150.0   5.3 224.0   0.00 22.61 26.72
mod3.1 2395.4   3.4 469.4   0.00 22.94 31.09

n.trials <- 1e4
age.seq <- seq(from = -2, to = 3.5, length.out = 58)
prediction.data <- data.frame(age = age.seq)

computeMu <- function(model, data, n.trials) {

    mu <- link(fit = model, data = data, n = n.trials)
    return(mu)

}

computeMuMean <- function(mu) {

    mu.mean <- apply(X = mu, MARGIN = 2, FUN = mean)
    return(mu.mean)

}

computeMuHPDI <- function(mu) {

    mu.HPDI <- apply(X = mu, MARGIN = 2, FUN = HPDI, prob = 0.97)
    return(mu.HPDI)

}

simulateHeights <- function(model, prediction.data) {

    simulated.heights <- sim(fit = model, data = prediction.data)
    return(simulated.heights)

}

plotResults <- function(model, prediction.data, original.data, n.trials) {

    mu <- computeMu(model, prediction.data, n.trials)
    mu.mean <- computeMuMean(mu)
    mu.HPDI <- computeMuHPDI(mu)
    simulated.heights <- simulateHeights(model = model, prediction.data = prediction.data)
    simulated.heights.HPDI <- apply(X = simulated.heights, MARGIN = 2, FUN = HPDI)
    plot(height ~ age, data = original.data, col = "steelblue", pch = 16)
    lines(x = prediction.data$age, y = mu.mean, lty = 2)
    lines(x = prediction.data$age, y = mu.HPDI[1, ], lty = 2)
    lines(x = prediction.data$age, y = mu.HPDI[2, ], lty = 2)
    shade(object = simulated.heights.HPDI, lim = prediction.data$age)

}

plotResults(
    model = mod3.1, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.2, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.3, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.4, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.5, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

plotResults(
    model = mod3.6, prediction.data = prediction.data,
    original.data = d1, n.trials = n.trials
    )

h.ensemble <- ensemble(mod3.4, mod3.5, mod3.6, data = list(age = age.seq) )
mu.mean <- apply(h.ensemble$link, 2, mean)
mu.ci <- t(apply(h.ensemble$link, 2, HPDI) )
height.ci <- t(apply(h.ensemble$sim, 2, HPDI) )

ggplot(data = d1, aes(x = as.numeric(age), y = height) ) +
    geom_point(size = 2) +
    geom_line(data = data.frame(age = age.seq, height = mu.mean) ) +
    geom_ribbon(
        data = data.frame(mu.ci), inherit.aes = FALSE,
        aes(x = age.seq, ymin = mu.ci[, 1], ymax = mu.ci[, 2]), alpha = 0.2) +
    geom_ribbon(
        data = data.frame(height.ci), inherit.aes = FALSE,
        aes(x = age.seq, ymin = height.ci[, 1], ymax = height.ci[, 2]), alpha = 0.1) +
    theme_bw(base_size = 20) + xlab("age") + xlim(-2, 3) + ylim(60, 190)

# model 1
coefs <- coef(mod3.1)
mu <- coefs["alpha"] + coefs["beta.1"] * d2$age
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.1 <- -2 * log.likelihood

# model 2
coefs <- coef(mod3.2)
mu <- coefs["alpha"] + coefs["beta.1"] * d2$age + coefs["beta.2"] * (d2$age)^2
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.2 <- -2 * log.likelihood

# model 3
coefs <- coef(mod3.3)
mu <- coefs["alpha"] + coefs["beta.1"] * d2$age + coefs["beta.2"] * (d2$age)^2 + coefs["beta.3"] * (d2$age)^3
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.3 <- -2 * log.likelihood

# model 4
coefs <- coef(mod3.4)
mu <- coefs["alpha"] + coefs["beta.1"]*d2$age + coefs["beta.2"]*(d2$age)^2 + coefs["beta.3"]*(d2$age)^3 + coefs["beta.4"]*(d2$age)^4
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.4 <- -2 * log.likelihood

# model 5
coefs <- coef(mod3.5)
mu <- coefs["alpha"] + coefs["beta.1"]*d2$age + coefs["beta.2"]*(d2$age)^2 + coefs["beta.3"]*(d2$age)^3 + coefs["beta.4"]*(d2$age)^4 + coefs["beta.5"]*(d2$age)^5
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.5 <- -2 * log.likelihood

# model 6
coefs <- coef(mod3.6)
mu <- coefs["alpha"] + coefs["beta.1"]*d2$age + coefs["beta.2"]*(d2$age)^2 + coefs["beta.3"]*(d2$age)^3 + coefs["beta.4"]*(d2$age)^4 + coefs["beta.5"]*(d2$age)^5 + coefs["beta.6"]*(d2$age)^6
log.likelihood <- sum(dnorm(x = d2$height, mean = mu, sd = coefs["sigma"], log = TRUE) )
dev.mod3.6 <- -2 * log.likelihood

deviances <- c(dev.mod3.1,dev.mod3.2,dev.mod3.3,dev.mod3.4,dev.mod3.5,dev.mod3.6)
comparison <- compare(mod3.1,mod3.2,mod3.3,mod3.4,mod3.5,mod3.6)
comparison <- as.data.frame(comparison@output)
comparison <- comparison[order(rownames(comparison) ), ]
waics <- comparison$WAIC

data.frame(deviance = deviances, waic = waics) %>%
    gather(type, value) %>%
    mutate(x = rep(1:6, 2) ) %>%
    ggplot(aes(x = x, y = value, colour = type) ) +
    scale_colour_grey() +
    geom_point(size = 2) +
    scale_x_continuous(breaks = 1:6) +
    theme_bw(base_size = 20) + xlab("model") + ylab("Déviance/WAIC")